
module predata

use params
use normcdf

implicit none

double precision wlthdat(2688,89)
double precision HHIDdat(2688)
double precision Mstatdat(2688,6)
double precision agedat(2688,6)
double precision asstdat(2688,6)

double precision num(nosims)

double precision, parameter :: PI(3)=(/0.17,0.5,0.83/) !median value for each of the three terciles (ie 1-33, 34-66, 67-99)

integer, parameter :: bornage = 66  !First year
integer, parameter :: dieage = 100  !Final possible year

integer, private :: i,h,w !counters
integer, private :: c
integer, parameter :: ashft=1
integer totobs
integer :: outcsv !output to csv files

contains

!------------------------------------------------------------------------------------------------------------------
subroutine readin

implicit none

!for the coef matrices:
! 	column 1: age
! 	column 2: age shifter
! 	column 3: health shifter
! 	column 4: male shifter
! 	column 5: permanent income shifter (multiply by percentile)
! 	column 6: permanent income squared (multiiply by percentile squared)

double precision hlthcoef(nh*(dieage-bornage)/2,2+4*(nh-1))
double precision hlthprob(nh*(dieage-bornage)/2,nwp,nh) !just by health and income tercile (assume female)

double precision deathcoef(nh*(dieage-bornage)/2,6)
double precision deathprob(nh*(dieage-bornage)/2,nwp) !just by health and income tercile (assume female)

!**********************************************************************************
!HEALTH PROBABILITIES (PROBABILITY OF BEING IN BAD HEALTH)
open(1,file='input/healthtransitions.out',status='old')
do i=1,nh*(dieage-bornage)/2
	read (1,*) hlthcoef(i,:)
enddo
close(1)
!1=healthy, 2=semi-unhealthy, 3=unhealthy
do i=1,(dieage-bornage)/2 !over age
	do h=1,nh !over health states yesterday
		do w=1,nwp
			c=(i-1)*nh+h
			!               ageshift            PIshift		    femaleshift		healthshift
			hlthprob(c,w,2)=hlthcoef(c,3)+PI(w)*hlthcoef(c,4)+hlthcoef(c,5)+hlthcoef(c,6)
			hlthprob(c,w,3)=hlthcoef(c,7)+PI(w)*hlthcoef(c,8)+hlthcoef(c,9)+hlthcoef(c,10)
			!multinomial logit:
			hlthprob(c,w,1)=							   1.0/(1.0+exp(hlthprob(c,w,2))+exp(hlthprob(c,w,3)))
			hlthprob(c,w,2)=exp(hlthprob(c,w,2))/(1.0+exp(hlthprob(c,w,2))+exp(hlthprob(c,w,3)))
			hlthprob(c,w,3)=exp(hlthprob(c,w,3))/(1.0+exp(hlthprob(c,w,2))+exp(hlthprob(c,w,3)))
			!convert to one-year transitions:
			!hlthprob(c,3)=sqrt(hlthprob(c,3))
			!hlthprob(c,4)=sqrt(hlthprob(c,4))
			!hlthprob(c,5)=sqrt(hlthprob(c,5))
		enddo
	enddo
enddo

!reshape this matrix into the matrix used in main program
do i=1,nt !age
	do h=1,nh !health yesterday
		do w=1,nwp
			if (annualperiods.eq.1) then
				if (MOD(i,2).eq.0) then 
					c=((i/2)-1)*nh+h
					hprob(1,h,w,i)=hlthprob(c,w,1)
					hprob(2,h,w,i)=hlthprob(c,w,2)
					hprob(3,h,w,i)=hlthprob(c,w,3)
				else ! assume same health for sure on off years
					hprob(:,h,w,i)=0.0
					hprob(h,h,w,i)=1.0 !override
				endif
			else
				c=(i-1)*nh+h
				hprob(1,h,w,i)=hlthprob(c,w,1)
				hprob(2,h,w,i)=hlthprob(c,w,2)
				hprob(3,h,w,i)=hlthprob(c,w,3)
	! 			print*,'health prob today given h,i,w',h,i,w,hprob(:,h,w,i)
			endif
		enddo
	enddo
enddo
!**********************************************************************************

!**********************************************************************************
!SURVIVAL PROBABILITIES
open(1,file='input/deathtransitions.out',status='old')
do i=1,nh*(dieage-bornage)/2
	read (1,*) deathcoef(i,:)
	!print*, deathcoef(i,:)
enddo
close(1)
!1=healthy, 2=semi-unhealthy, 3=unhealthy
do i=1,(dieage-bornage)/2 !over age
	do h=1,nh !over health states yesterday
		do w=1,nwp
			c=(i-1)*nh+h
			!               age                  PIshift        femaleshift	   healthshift
			deathprob(c,w)=deathcoef(c,3)+PI(w)*deathcoef(c,4)+deathcoef(c,5)+deathcoef(c,6)
			!logit:
			deathprob(c,w)=exp(deathprob(c,w))/(1+exp(deathprob(c,w)))
			!print*,deathprob(c,:)
		enddo
	enddo
enddo

! !reshape this matrix into the matrix used in main program
!do i=1,dieage-bornage !age
do i=1,nt !age
	do h=1,nh !health yesterday
		do w=1,nwp
			if (annualperiods.eq.1) then
				if (MOD(i,2).eq.0) then 
					c=((i/2)-1)*nh+h
					survprob(h,w,i)=1.0-deathprob(c,w)
				else
					survprob(h,w,i)=1.0 ! assume survive for sure on off years
				endif
			else
				c=(i-1)*nh+h
	! 			if (i.eq.nt) then
	! 				survprob(h,w,i)=0.0
	! 			else
	! 				survprob(h,w,i)=1.0-deathprob(c,w)
	! 			endif
				survprob(h,w,i)=1.0-deathprob(c,w)
			endif
			!print*,'surv prob today given h,i,w',h,i,w,survprob(h,w,i)
		enddo
	enddo
enddo
!**********************************************************************************

!**********************************************************************************
!LOAD HRS PROFILES

if (annualperiods.eq.1) then
	open(1,file='input/input_profiles_cohortannual.out',status='old')
else 
	open(1,file='input/input_profiles_cohort.out',status='old')
endif
do i=1,nHRS
read(1,*) HRSprof_cohort(i),HRSprof_age(i),HRSprof_passets(i),HRSprof_pinc(i),HRSprof_ltci(i), &
		  HRSprof_kinc(i),HRSprof_kassets(i),HRSprof_health(i,:)
enddo

!**********************************************************************************

end subroutine readin
!------------------------------------------------------------------------------------------------------------------

!------------------------------------------------------------------------------------------------------------------
subroutine readin_moments

	implicit none

	open(1,file='input/parent_asset_cohort1_moments.out',status='old')
	do i=1,7
		read (1,*) passet_c1_p25_data(i,:),passet_c1_p50_data(i,:),passet_c1_p75_data(i,:)
	enddo
	close(1)
	open(1,file='input/parent_asset_cohort2_moments.out',status='old')
	do i=1,7
		read (1,*) passet_c2_p25_data(i,:),passet_c2_p50_data(i,:),passet_c2_p75_data(i,:)
	enddo
	close(1)
	open(1,file='input/parent_asset_cohort3_moments.out',status='old')
	do i=1,7
		read (1,*) passet_c3_p25_data(i,:),passet_c3_p50_data(i,:),passet_c3_p75_data(i,:)
	enddo
	close(1)
	open(1,file='input/parent_asset_cohort4_moments.out',status='old')
	do i=1,7
		read (1,*) passet_c4_p25_data(i,:),passet_c4_p50_data(i,:),passet_c4_p75_data(i,:)
	enddo
	close(1)
	open(1,file='input/parent_asset_cohort5_moments.out',status='old')
	do i=1,7
		read (1,*) passet_c5_p25_data(i,:),passet_c5_p50_data(i,:),passet_c5_p75_data(i,:)
	enddo
	close(1)	

	
	open(1,file='input/child_asset_cohort1_moments.out',status='old')
	do i=1,7
		read (1,*) kasset_c1_p25_data(i,:),kasset_c1_p50_data(i,:),kasset_c1_p75_data(i,:)
	enddo
	close(1)
	open(1,file='input/child_asset_cohort2_moments.out',status='old')
	do i=1,7
		read (1,*) kasset_c2_p25_data(i,:),kasset_c2_p50_data(i,:),kasset_c2_p75_data(i,:)
	enddo
	close(1)
	open(1,file='input/child_asset_cohort3_moments.out',status='old')
	do i=1,7
		read (1,*) kasset_c3_p25_data(i,:),kasset_c3_p50_data(i,:),kasset_c3_p75_data(i,:)
	enddo
	close(1)
	open(1,file='input/child_asset_cohort4_moments.out',status='old')
	do i=1,7
		read (1,*) kasset_c4_p25_data(i,:),kasset_c4_p50_data(i,:),kasset_c4_p75_data(i,:)
	enddo
	close(1)
	open(1,file='input/child_asset_cohort5_moments.out',status='old')
	do i=1,7
		read (1,*) kasset_c5_p25_data(i,:),kasset_c5_p50_data(i,:),kasset_c5_p75_data(i,:)
	enddo
	close(1)

end subroutine readin_moments
!------------------------------------------------------------------------------------------------------------------

!------------------------------------------------------------------------------------------------------------------
subroutine readin_varcov 

	implicit none

	open(1,file='input/variance_covariance_matrix.out',status='old')
	do i=1,94 !Nmoments
		!read (1,*) varcov_matrix(i,1:Nmoments)
		read (1,*) varcov_matrix(i,1:94)
	enddo
	close(1)

end subroutine readin_varcov
!------------------------------------------------------------------------------------------------------------------

subroutine life_expectancy_table

	implicit none

	integer index, node, t
	integer HRSperson
	integer PI
	integer simhealth(nosims,nt)
	integer startage(nosims)

	do index=1,nosims

		HRSperson=ceiling(dble(nHRS)*randprofile(index))
		PI=HRSprof_pinc(HRSperson)
		startage(index)=HRSprof_age(HRSperson)
		simhealth(index,1)=HRSprof_health(HRSperson,1)

		do t=2,nt
			node=(index-1)*(nt) + t
			
			!first check if dead already
			if (simhealth(index,t-1)==0) then
				simhealth(index,t)=0
			!if not dead already check if dead now
			elseif (randsurv(node)>=survprob(simhealth(index,t-1),PI,t)) then
				simhealth(index,t)=0
			!if not dead find health
			elseif (randhlth(node)<=hprob(1,simhealth(index,t-1),PI,t)) then
				simhealth(index,t)=1
			elseif (randhlth(node)<=(hprob(1,simhealth(index,t-1),PI,t) + &
															 hprob(2,simhealth(index,t-1),PI,t))) then
				simhealth(index,t)=2
			else
				simhealth(index,t)=3
			endif

		enddo !over time

	enddo !over index of simulated people

	!LIFE EXPECTANCY/HEALTH TABLE
	if (pid==0) then
	  open(1,file='output/simhealth_fromestimated.csv',action="write",status="replace")
		do index=1,nosims
		    write(1,"(100000(I3,',',:))") startage(index),simhealth(index,:)
		enddo
	endif

end subroutine life_expectancy_table

end module predata